Regresión lineal multiple y multicolinealidad
Juan Jose Rivera Londoño
29/8/2021
Basede datos.
A continuación se leerá la base de datos sobre la eficacia en el control de infecciones hospitalarias.
data = read.table("D:/SEMESTRE_1_2021/ESTADISTICA_2/INFORME_2/datos2.txt", header = TRUE, sep = " ", dec = ".")1. Estime un modelo de regresión lineal múltiple que explique el Riesgo de Infección en términos de todas las variables predictoras. Analice la significancia de la regresión y de los parámetros individuales. Interprete los parámetros estimados. Calcule e interprete el coeficiente de determinación múltiple R2. Comente los resultados.
A continuación se estimará un modelo con todas las variables predictoras.
##
## Call:
## lm(formula = Y ~ ., data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.03944 -0.80043 -0.00266 0.60450 2.23292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5986009 1.5159559 -0.395 0.694365
## X1 0.2106683 0.0785765 2.681 0.009501 **
## X2 0.0197512 0.0277108 0.713 0.478803
## X3 0.0470925 0.0132888 3.544 0.000779 ***
## X4 0.0105604 0.0073166 1.443 0.154213
## X5 0.0008996 0.0007379 1.219 0.227679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.035 on 59 degrees of freedom
## Multiple R-squared: 0.4468, Adjusted R-squared: 0.3999
## F-statistic: 9.53 on 5 and 59 DF, p-value: 1.058e-06La prueba F para la regresión nos da un p-value igual a \(1.05 \times 10^-6\), con un nivel de significancia del 5% podemos decir que el modelo obtenido con la combinación de estas variables predictoras es significativo, esto no significa que todas las variables sean significativas, de hecho, solo las variables x1 (Duración de la estadía) y x3 (Número de camas) lo son, los demás coeficientes de regresión pueden ser cero dependiendo de una muestra diferente, lo cual implica que estadísticamente estas variables no sean significativas y no aporten información predictiva para Y (Riesgo de infección). Interpretaremos los coeficientes \(\beta_1\) y \(\beta_3\), el primero significa que en promedio el riesgo de infección incrementa en un 21% por un aumento en un día de estadía cuando las demás variables explicativas permanecen constantes, y el segundo coeficiente representa que el riesgo de infección aumenta en un 4% en promedio por un aumento de 1 cama en el hospital, también con las demás variables constantes. El coeficiente \(R^2\) del modelo es de \(0.44\) esto significa que un 44% de la variabilidad total en la variable riesgo de infección es explicada por el modelo de regresión. A mi parecer la selección de variables predictoras no aporta mucha información para predecir el riesgo de infección.
2. Use la tabla de todas las regresiones posibles, para probar la significancia simultánea del subconjunto de tres variables con los valores p mayores del punto anterior. Según el resultado de la prueba es posible descartar del modelo las variables del subconjunto?.
El subconjunto de variables que tienen los mayores p-values son X2, X4, X5 con p-values iguales a 0.47, 0.15, 0.22 respectivamente. las hipótesis para este tipo de prueba es:
\[H_o: \beta_2=\beta4=\beta_5=0\] \[H_1: \text{alguno entre } \beta_2,\beta4,\beta_5 \neq 0\]
Para formular esta prueba de hipótesis se requieren 2 modelos, el primero es el modelo de regresión lineal con todas las variables predictoras y el segundo es el modelo reducido, en el cual se eliminan del modelo el subconjunto de variables de la hipótesis nula. El estadístico de prueba será entonces:
\[F=\frac{[SSE(X1,X3)-SEE(X1,X2,X3,X4,X5)]/[(n-3)-(n-6)]}{MSE(X1,X2,X3,X4,X5)}\]
O equivalentemente:
\[F=\frac{[SSR(X1 ,X2, X3, X4, X5)-SSR(X1,X3)]/3}{MSE(X1,X2,X3,X4,X5)}\]
Usando el código desarrollado por el profesor para obtener la tabla de todas las regresiones posibles, tenemos que el SSE del modelo completo es \(63.203\) y que el SSE para el modelo reducido es \(67.257\), también tenemos que el MSE del modelo completo es \(63.203/(n-6)\), donde \(n=65\), También se planteara un procedimiento alternativo pero equivalente sin usar la tabla de todas las regresiones posibles y calculando el SSR en vez del SSE, tenemos usando el SSE:
\[F=\frac{[67.257-63.203]/[(65-3)-(65-6)]}{63.203/(65-6)}=1.2614\]
Usando el código equivalente o el SSR encontramos que el F0 es el mismo.
mod2 = lm(Y~X1+X3,data=data)
myAllRegTable(mod1,MSE = FALSE)## k R_sq adj_R_sq SSE Cp Variables_in_model
## 1 1 0.249 0.237 85.813 19.106 X3
## 2 1 0.244 0.232 86.376 19.632 X1
## 3 1 0.182 0.169 93.468 26.252 X4
## 4 1 0.068 0.053 106.508 38.425 X5
## 5 1 0.001 -0.014 114.090 45.503 X2
## 6 2 0.411 0.392 67.257 3.784 X1 X3
## 7 2 0.319 0.297 77.804 13.629 X3 X4
## 8 2 0.311 0.289 78.739 14.503 X1 X4
## 9 2 0.293 0.270 80.764 16.392 X3 X5
## 10 2 0.276 0.253 82.667 18.169 X2 X3
## 11 2 0.258 0.235 84.721 20.087 X1 X5
## 12 2 0.248 0.224 85.877 21.166 X1 X2
## 13 2 0.234 0.209 87.519 22.698 X4 X5
## 14 2 0.182 0.156 93.416 28.203 X2 X4
## 15 2 0.071 0.041 106.102 40.046 X2 X5
## 16 3 0.430 0.402 65.108 3.778 X1 X3 X4
## 17 3 0.422 0.393 66.088 4.693 X1 X3 X5
## 18 3 0.415 0.386 66.850 5.404 X1 X2 X3
## 19 3 0.359 0.327 73.279 11.405 X3 X4 X5
## 20 3 0.336 0.303 75.868 13.822 X2 X3 X4
## 21 3 0.328 0.295 76.791 14.684 X1 X4 X5
## 22 3 0.325 0.292 77.090 14.963 X2 X3 X5
## 23 3 0.314 0.280 78.400 16.186 X1 X2 X4
## 24 3 0.261 0.224 84.459 21.842 X1 X2 X5
## 25 3 0.236 0.198 87.321 24.514 X2 X4 X5
## 26 4 0.442 0.405 63.748 4.508 X1 X3 X4 X5
## 27 4 0.433 0.395 64.795 5.486 X1 X2 X3 X4
## 28 4 0.427 0.389 65.435 6.083 X1 X2 X3 X5
## 29 4 0.379 0.338 70.904 11.188 X2 X3 X4 X5
## 30 4 0.329 0.284 76.656 16.558 X1 X2 X4 X5
## 31 5 0.447 0.400 63.203 6.000 X1 X2 X3 X4 X5n=dim(data)[1]
dfn=6-3
dfd=n-6
f0=((67.257-63.203)/((65-3)-(65-6)))/(63.203/(65-6))
SRR1 = sum((fitted.values(mod1)-mean(data$Y))**2)
SRR2 = sum((fitted.values(mod2)-mean(data$Y))**2)
SSE1 = sum((fitted.values(mod1)-data$Y)**2)
MSE1 = SSE1/(n-6)
F0 = ((SRR1-SRR2)/3)/(MSE1)
## Cuantil de la distribucion f
qf(1-0.05,dfn,dfd)## [1] 2.760767## Estadistico de prueba con SSE
f0## [1] 1.26147## Estadistico de prueba con SSR
F0## [1] 1.261376Usando un nivel de significancia del 5%, Como el valor critico o cuantil \(F_{\alpha,3,59}=2.76\) tenemos:
\[F_0<F_{\alpha,3,59}\ ; \ 1.26<2.76 \]
Por lo tanto no podemos rechazar la hipótesis nula y el subconjunto de variables X2, X4, X5 no son estadisticamente significativas, por lo tanto si es posible descartarlas del modelo.
3. Plantee una pregunta donde su solución implique el uso exclusivo de una prueba de hipótesis lineal general de la forma \(H_0 : L{\times} \beta = 0\) (solo se puede usar este procedimiento y no SSextra), donde especifique claramente la matriz L, el modelo reducido y la expresión para el estadístico de prueba.
Evaluaremos simultáneamente las hipótesis: \[\beta_2 = \beta_4\] \[\beta_1 = \beta_3\] estas pueden ser reescritas como: \[\beta_2 - \beta_4 = 0\] \[\beta_1 - \beta_3 = 0\] O en forma matricial si consideramos todos los parámetros del modelo: \[ \left[\begin{array}{cc} 0 & 0 & 1 & 0 & -1 & 0\\ 0 & 1 & 0 & -1 & 0 & 0 \end{array}\right] \times \left[\begin{array}{cc} \beta_0\\ \beta_1\\ \beta_2\\ \beta_3\\ \beta_4\\ \beta_5 \end{array}\right] = \left[\begin{array}{cc} 0\\ 0 \end{array}\right]\] Donde la matriz L es igual a: \[ L = \left[\begin{array}{cc} 0 & 0 & 1 & 0 & -1 & 0\\ 0 & 1 & 0 & -1 & 0 & 0 \end{array}\right] \] Para probar la hipótesis dada anteriormente se plantean 2 modelos, uno en el cual la hipótesis nula se cumple, o modelo reducido y el otro es el modelo completo, para el modelo reducido tenemos: \[ Y = \beta_0+\beta_1X_{1,3}+\beta_2X_{2,4}+\beta_5X_5 \] Donde: \[ X_{1,3}=X_1+X_3 \qquad \text{y} \qquad X_{2,4}=X_2+X_4 \]
Ahora teniendo los 2 modelos de regresión el estadístico de prueba es el siguiente:
\[ F_0 = \frac{SSR(X1,X2, X3, X4, X5)-SSR(X13, X24, X5)/(5-3)}{MSE(X1, X2, X3, X4, X5)} \] Donde los grados de libertad del numerador son \(5-3=2\) y los grados de libertad del denominador son \(n-6=59\).
y=data$Y
x13 =data$X1+data$X3
x24 =data$X2+data$X4
x5=data$X5
dfnl=2
dfdl=n-6
mod3 = lm(y~x13+x24+x5)
SRR3 = sum((fitted.values(mod3)-mean(data$Y))**2)
F0l = ((SRR1-SRR3)/dfnl)/(MSE1)
qf2 = qf(1-0.05,df1 = dfnl,df2=dfdl)
## Estadístico de prueba
F0l## [1] 2.345988##Valor crítico de la distribucion F
qf2## [1] 3.153123El estadístico de prueba para la hipótesis planteada es \(F_0 = 2.34\), Con un nivel de significáncia del 5% el valor critico de prueba de la distribución F es \(F_{\alpha,2,59} = 3.15\) por lo tanto.
\[F_0<F_{\alpha,2,59} \qquad ; \qquad 2.34<3.15\]
Por lo tanto no podemos rechazar la hipótesis nula.
4. Realice una validación de los supuestos en los errores y examine si hay valores atípicos, de balanceo e influenciales. Qué puede decir acerca de la validez de éste modelo?. Argumente.
Para validar los supuestos de los errores tenemos que revisar que los residuales:
- Tengan varianza constante
- Se distribuyan aproximadamente normal
- Sean independientes
- Tengan media cero
Pero ahora lo haremos con los residuales estudentizados ya que deberíamos escalarlos para poder determinar verdaderamente cuando un residual es “grande”.
Varianza constante
Este supuesto se puede validar con un gráfico de residuales estudentizados del modelo vs valores ajustados de la variable respuesta.
res_vs_fit = as.data.frame(cbind(fitted.values(mod1),studres(mod1)))
plt1 = ggplot(data=res_vs_fit,aes(V1,V2)) +
geom_point()+theme_bw()+
labs(x="Valores ajustados",y = "Residuales estudentizados")
ggplotly(plt1)Como podemos ver en la gráfica la distribución de los residuales es homogénea en todo el espacio, hay una aglomeración de puntos en la parte izquierda del gráfico pero esto es por la cantidad de predicciones que caen en ese rango mas no porque los residuales varían conforme aumenta o disminuyan los valores ajustados, podemos decir que la varianza es constante.
Distribución normal en los errores.
Para validar este supuesto podemos realizar y analizar el gráfico de probabilidad normal de los residuales.
qqnorm(studres(mod1),
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Normality test")
qqline(studres(mod1))
shapiro.test(studres(mod1))##
## Shapiro-Wilk normality test
##
## data: studres(mod1)
## W = 0.98663, p-value = 0.7094Como podemos observar en el gráfico los puntos mas alejados del centro de la gráfica (en especial los puntos a la izquierda) se alejan de la recta de normalidad, cabe destacar que no son muchos los puntos que se alejan de esta recta y ademas la prueba de Shapiro-Wilk nos da un p-value igual a 0.7, por esta razón concluimos que los errores si siguen una distribución normal. Aunque concluimos que los residuales si distribuyen normal los puntos a la izquierda del gráfico pueden corregirse a menudo con una transformación de los valores de Y , esto logra estabilizar la varianza y una mayor aproximación a la normalidad.
Independencia y media cero.
Como no tenemos un registro temporal o una secuencia en la cual fueron tomados los datos no podemos concluir nada acerca de la independencia de los errores o no tiene sentido en este caso concluir acerca de este hecho.
Por la manera de conformar el modelo lineal, ya sea por mínimos cuadrados o por máxima verosimilitud (son equivalentes en la determinación de los coeficientes y en esta propiedad en especifico), los residuales del modelo siempre tienen media cero.
Observaciones atipicas
Para analizar si en el modelo se encuentran observaciones atípicas, se observan los residuales estudentizados del modelo, estos se calculan como:
res_stud = studres(mod1)
res_stan = rstandard(mod1)
res_stud[abs(res_stud)>3]## named numeric(0)res_stan[abs(res_stan)>3]## named numeric(0)Como podemos ver del código anterior, no hay ninguna observación para las cuales se observen que \(|ri|\) o \(|di|\) sea mayor a 3, en este caso no tenemos observaciones atípicas.
Puntos de balanceo.
Los puntos de balanceo son detectados mediante el análisis de los elementos de la diagonal principal de la matriz \(H\), los \(hii\), se decide si una observación es un punto de balanceo si \(hii>\frac{2p}{n}\), donde p es el numero de parámetros del modelo y \(\bar{hii}=p/n\) es la media de los \(hii\), si \(2p=n > 1\) este criterio no funciona pues los hii siempre son menores a 1.
x0 = rep(1,dim(data)[1])
X = matrix(c(x0,as.matrix(data[,2:6])),ncol=6)
colnames(X) = c('x0','x1','x2','x3','x4','x5')
betas = solve(t(X) %*% X)%*%(t(X)%*%as.matrix(data$Y))
H = X%*%solve(t(X)%*%X)%*%t(X)
Hii = as.matrix(diag(H))
Hb = dim(betas)[1]/n
## Media de hii
Hb## [1] 0.09230769## valores atipicos
atip = as.matrix(Hii[abs(Hii)>2*Hb])
row.names(atip) = which(Hii>=2*Hb)
atip## [,1]
## 22 0.2006660
## 23 0.2242422
## 27 0.4471487
## 30 0.2013972
## 45 0.2348574
## 52 0.3432044como podemos ver \(\bar{hii}=0.09\), por lo que esta prueba es valida y también podemos ver que las observaciones 22, 23, 27, 30, 45 y 52 son mayores a \(2\times Hb = 0.18\) por lo tanto son puntos de balanceo.
Puntos influyentes.
Para probar si una observación es un punto influyente utilizamos distintas medidas una es la distancia de Cook y la otra es la distancia DFFITS.el punto sera influyente si la distancia de Cook \(Di>1\) o si la distancia \(|DFFITS_i|>2 \sqrt{\frac{p}{n}}\).
p = dim(betas)[1]
cd = cooks.distance(mod1)
DFFI = dffits(mod1)
### Puntos influyentes
cd[cd>1]## named numeric(0)DFFI[abs(DFFI)>2*sqrt(p/n)]## 30 41 45 52
## -0.6585742 0.7429419 0.7849038 -1.4021459### DFBETAS para cada variableComo podemos ver bajo el criterio de la distancia de Cook no se detecto ningún punto influyente pero bajo el criterio de la distancia DFFITS se detectó que los puntos 30, 41, 45 y 52 son puntos influyentes.
Concluimos finalmente que el modelo si tiene validez en cuanto a los supuestos de los errores ya que estos se cumplen, aun así hay varios puntos que son atípicos e influenciales en el modelo, el \(R^2\) del modelo es bajo a mi parecer, los puntos de balanceo tienen una influencia directa en esta estadística de desempeño del modelo, un chequeo de estos puntos es importante ya que nos indican puntos que pueden ser eliminados del modelo llegado el caso de que sean errados (errores de medición etc), los puntos de balanceo “halan” y modifican la regresión en su dirección, se sugiere que se chequeen estos puntos para saber si verdaderamente representan la realidad o si hay errores errores de medición en los datos.
5. Verificar la presencia de multicolinealidad usando graficos y/o indicadores apropiados.
Se verificará la presencia de multicolinealidad por medio de un gráfico de la matriz de correlación entre las diferentes variables predictoras y haciendo uso de indicadores como los factores de inflación de la varianza, indices de condición y proporciones de descomposición de varianza.
Gráfico de correlación.
corr_mat = cor(data[,2:6])
corrplot(corr_mat, type="upper", order="hclust", tl.col="black", tl.srt=45)
pairs(data,lower.panel =myPanel.cor,upper.panel = panel.smooth,diag.panel=myPanel.box,
labels=names(data))
Como se observa en la anterior gráfica, la correlación entre variables predictoras no es muy fuerte, el valor critico para determinar indicios de correlación entre las predictoras es de 0.5 en valor absoluto, como podemos ver el valor máximo es entre las variables es de 0.38 entre x1 y x4, por lo tanto esta prueba no nos da indicios de correlación, y además no es concluyente.
Factor de inflación de varianza.
Utilizaremos la función VIF() para calcular la inflación de la varianza, si el VIF es mayor a 5 no tenemos problemas de multicolinealidad.
##Según Vif no se encuentran efectos de multicolinealidad
vf = vif(mod1)
vf## X1 X2 X3 X4 X5
## 1.375666 1.176970 1.270949 1.302309 1.124675## Como podemos ver el valor para cada variable es menor a 5 por lo que no detectamos problemas de multicolinealidad con esta prueba.
Indice de condición y proporciones de varianza.
Calcularemos los indices de condición, si el indice de condición es menor a 10, no tenemos multicolinealidad, si el indice de condición esta entre 10 y 31.62 tenemos multicolinealidad moderada, tenemos un indice de condición mayor a 31.2 tenemos multicolinealidad severa, también analizaremos las proporciones de la varianza \(\pi_{i,j}\) si este valor es mayor a 0.5 para 2 o mas variables asociadas a un mismo valor propio podemos decir que existen problemas de multicolinealidad entre las respectivas variables.
x1_stan = (data$X1-mean(data$X1))/sqrt(sum((data$X1-mean(data$X1))^2)/(n-1))
x2_stan = (data$X2-mean(data$X2))/sqrt(sum((data$X2-mean(data$X2))^2)/(n-1))
x3_stan = (data$X3-mean(data$X3))/sqrt(sum((data$X3-mean(data$X3))^2)/(n-1))
x4_stan = (data$X4-mean(data$X4))/sqrt(sum((data$X4-mean(data$X4))^2)/(n-1))
x5_stan = (data$X5-mean(data$X5))/sqrt(sum((data$X5-mean(data$X5))^2)/(n-1))
Y_stan = (data$Y-mean(data$Y))/sqrt(sum((data$Y-mean(data$Y))^2)/(n-1))
std_mod = lm(Y_stan~x1_stan+x2_stan+x3_stan+x4_stan+x5_stan)
summary(std_mod)##
## Call:
## lm(formula = Y_stan ~ x1_stan + x2_stan + x3_stan + x4_stan +
## x5_stan)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.52642 -0.59908 -0.00199 0.45244 1.67124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.788e-17 9.608e-02 0.000 1.000000
## x1_stan 3.045e-01 1.136e-01 2.681 0.009501 **
## x2_stan 7.488e-02 1.051e-01 0.713 0.478803
## x3_stan 3.869e-01 1.092e-01 3.544 0.000779 ***
## x4_stan 1.595e-01 1.105e-01 1.443 0.154213
## x5_stan 1.252e-01 1.027e-01 1.219 0.227679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7747 on 59 degrees of freedom
## Multiple R-squared: 0.4468, Adjusted R-squared: 0.3999
## F-statistic: 9.53 on 5 and 59 DF, p-value: 1.058e-06eigprop(mod1, na.rm = TRUE, Inter = TRUE, prop = 0.5)##
## Call:
## eigprop(mod = mod1, na.rm = TRUE, Inter = TRUE, prop = 0.5)
##
## Eigenvalues CI (Intercept) X1 X2 X3 X4 X5
## 1 5.4018 1.0000 0.0002 0.0009 0.0002 0.0064 0.0014 0.0079
## 2 0.3026 4.2251 0.0002 0.0001 0.0002 0.1557 0.0034 0.7766
## 3 0.2365 4.7788 0.0024 0.0050 0.0038 0.6408 0.0043 0.1243
## 4 0.0343 12.5430 0.0223 0.0037 0.0285 0.1144 0.9209 0.0021
## 5 0.0208 16.1205 0.0399 0.9761 0.0293 0.0042 0.0654 0.0677
## 6 0.0039 37.1276 0.9349 0.0141 0.9378 0.0785 0.0047 0.0214
##
## ===============================
## Row 5==> X1, proportion 0.976147 >= 0.50
## Row 6==> X2, proportion 0.937849 >= 0.50
## Row 3==> X3, proportion 0.640772 >= 0.50
## Row 4==> X4, proportion 0.920851 >= 0.50
## Row 2==> X5, proportion 0.776606 >= 0.50eigprop(std_mod, na.rm = TRUE, Inter = TRUE, prop = 0.5)##
## Call:
## eigprop(mod = std_mod, na.rm = TRUE, Inter = TRUE, prop = 0.5)
##
## Eigenvalues CI (Intercept) x1_stan x2_stan x3_stan x4_stan x5_stan
## 1 1.7324 1.0000 0 0.1327 0.0002 0.1074 0.1421 0.0648
## 2 1.2397 1.1821 0 0.0844 0.4434 0.1317 0.0011 0.0000
## 3 1.0000 1.3162 1 0.0000 0.0000 0.0000 0.0000 0.0000
## 4 0.9875 1.3245 0 0.0096 0.0283 0.0566 0.1361 0.6375
## 5 0.5461 1.7810 0 0.0136 0.2065 0.6457 0.5726 0.0027
## 6 0.4943 1.8720 0 0.7597 0.3217 0.0586 0.1481 0.2950
##
## ===============================
## Row 6==> x1_stan, proportion 0.759749 >= 0.50
## Row 5==> x3_stan, proportion 0.645686 >= 0.50
## Row 5==> x4_stan, proportion 0.572581 >= 0.50
## Row 4==> x5_stan, proportion 0.637512 >= 0.50Según los indices de condición calculados, tenemos 2 problemas de multicolinealidad moderada y uno de alta multicolinealidad. Analizando los indices de descomposición de la varianza del modelo sin estandarizar no podemos determinar que variables tienen los problemas de multicolinealidad, pero si miramos la descomposición del modelo estandarizado podemos que las variables X3 Y X4 presentan presentan un problema de multicolinealidad moderado esto se puede ver en el valor propio numero 5 ya que el indice de condición de este valor propio en el modelo sin estandarizar es igual a 16.12, se debería analizar si eliminar una de las variables para eliminar los problemas asociados.